function output=discrete(n,p,Sup)
%Step 1: find LCM using p 
ps=sym(p); %ax1
[~, den]=numden(ps);
dendouble=double(den); %ax1
common=elcm(dendouble); %1x1

%Step 2: draw n numbers from Uniform {1,...,common}
Ind=randi(common,n,1); %nx1

%Step 3: compute numerator as p*common
num=p*common;  %1xa

%Step 4: relabel each draws
output=zeros(n,1);
for i=1:n
    if Ind(i)<=num(1)
        output(i)=Sup(1);
    elseif Ind(i)>=common-num(end)+1
        output(i)=Sup(end);
    end
    for j=2:size(p,1)-1
        cs=cumsum(num(1:j));
        if Ind(i)>=cs(end-1)+1 && Ind(i)<=cs(end)
            output(i)=Sup(j);
        end
    end
end

end
